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Static and dynamic properties of matter-wave solitons in dense Bose-Einstein conden- 
sates, where three-body interactions play a significant role, have been studied by a vari- 
ational approximation (VA) and numerical simulations. For experimentally relevant pa- 
rameters, matter-wave solitons may acquire a flat-top shape, which suggests employing 
a super-Gaussian trial function for VA. Comparison of the soliton profiles, predicted by 
VA and those found from numerical solution of the governing Gross-Pitaevskii equation 
shows good agreement, thereby validating the proposed approach. 

Keyviords: Bose-Einstein condensat; flat-top soliton; variational approach 

1. Introduction 

The existence and properties of solitons in Bose-Einstein condensates (BEC) has 
been the subject of considerable interest over the recent decade (for a review see 
articles ^ and books l^^) . All of the main types of matter- wave solitons, such as 
darkS brightEl, and gapEl, have been observed in the experiments. Dark solitons 
emerge in a BEC with repulsive interactions between atoms (nonlinearity is defo- 
cusing), while for the existence of bright solitons the interatomic interaction has 
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to be attractive (nonlinearity is focusing). Gap solitons develop in repulsive con- 
densates loaded in periodic potential of the optical lattice. While in homogeneous 
condensates solitons exist due to the balance between nonlinearity and dispersion, 
in a repulsive BEG subject to a periodic potential of the optical lattice, gap solitons 
come out of the interplay between nonlinearity and periodicity of the medium. 

Nowadays the conditions for the existence of these main types of solitons are 
usually created in magnetic and optical traps for cold quantum gases Meantime, 
the experimental techniques for trapping and manipulation of BEC's are progressing 
and novel conditions for the localized matter-waves are emerging. In this regard, 
development of atom chips for BEC's ^ and further advances in Feshbach reso- 
nance management of atomic interactions ^ can be mentioned as two important 
examples. Specifically, the atom chi p technology combined with techniques to sup- 
press three-body recombination 1^ allows to produce long-lived condensates with 
increased density, where the contribution of three-body scattering is dominant. On 
the other hand, even in condensates with normal density, inhibition of the two-body 
s-wave scattering length by a Feshbach resonance technique gives rise to enhanced 
role of the three-body effects. In these conditions the governing mean-field Gross- 
Pitaevskii equation (GPE) for the dynamics of BEG includes higher order nonlinear 
terms, alongside with the usual cubic term. As a consequence, localized matter- 
waves exhibit novel features both in their shape and dynamic behavior, compared 
to matter-wave solitons of the conventional GPE. 

In this work we focus on a specific class of bright solitons, so called "flat-top" 
solitons, which remain less explored in the context of BEG. Flat-top solitons in 
BEG's emerge when the repulsion between atoms, originating from three-atom col- 
lisions, prevails the attraction resulting from two-atom interactions. In terms of the 
governing GPE, this implies that the dcfocusing quintic nonlinearity is stronger than 
the focusing cubic nonlinearity. In the experiments such a situation is realized when 
the density of the condensate is high, for instance in BEG's on atom chips, or when 
the usual two-body interactions are suppressed by means of a Feshbach resonance 
technique. The properties BEG described by GPE with cubic-quintic nonlinearity 
and generic trap potential can be explored using the Lagrangian formalism. Our 
objective is the development of a Lagrangian based variational approximation (VA) 
for flat-top solitons using the super-Gaussian trial function. Static version of the 
VA provides stationary shape of the soliton, while the time-dependent version can 
be used for studying small amplitude oscillations around stationary states. 

The paper is structured as follows. In section II we formulate the model and 
present the governing equations. In sections III and IV, respectively, the static 
and dynamic versions of the VA have been developed. In concluding section V we 
formulate our main findings. 
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2. The model and main equations 

The dynamics of a quasi- ID Bose- Einstein condensate with two- and three-body 
interactions, confined to an exter nal tra p potential, is described by the following 
Gross-Pitaevskii equation (see e.g. I^^ l -*^^ !) 

where iplx^t) is the mean-field wave function of the condensate, m is the atomic 
mass, gi = Anh^as/m is the coefficient of cubic nonlinearity with as being the atomic 
s-wave scattering length, characterizing the strength of two-body interactions, 52 is 
the coefficient of quintic nonlinearity responsible for three-body interactions, aj_ = 
yjh/muj_i is the transverse harmonic oscillator length, TV is the number of atoms 
comprising the BEG, and V{x) — 'muj'^x'^/2 is the harmonic trap potential in the 
axial direction. Strengths of the trap potential in the axial and radial directions are 
given, respectively, through the trap frequencies uj^ and lu_i_. 

In deriving the Eq. ([1]) from the original 3D GPE it is assumed that the 
transverse confinement of the condensate is strong enough, so that its radial de- 
grees of freedom are fixed. In this condition characteristic energies of the trans- 
verse excitations huj± are much greater than the energy from the nonlinear term 
Eint — 4:Tr fi^ a sTio/m, where as,no,m are the atomic s-wave scattering length, peak 
density of the condensate and the atomic mass, respectively. Therefore, the dynamic 
evolution of the condensate is possible only in the axial direction, while in the ra- 
dial direction it remains in the ground state of the strong parabolic trap (transverse 
dynamics is frozen) . Alternatively, the criterion for the one-dimensionality of BEG 
can be expressed as ^ h/muj± >> ^ = {Annoas)^^^'^ , which can be understood as 
the radial harmonic oscillator length being much greater than the healing length 
(for ^''Rb: as = 5.6 nm, no = 10^"^ cm~^, m = 1.45 x lO'^^fcg, £, = 0.4/ito). 
The wave function is subject to the normalization condition, 

f +00 

\ip{x)fdx^l. (2) 



For convenience we reduce the Eq. ([T]) into dimensionless form by adopting new 
variables t — )■ uj±t, x — > x/a±, a = —giN / (2T:a\ hLO±), 
(3 = 52TVV(37r2ai fii^^), U{x) = -V{x)/{nw^), 4, ^ 

+ ^V-.. + Uix)iP + alVPV - m^i^ = 0. (3) 

It is instructive to estimate the experimentally realistic values for the main pa- 
rameters a and /3 in Eq. According to data of for ®^Rb condensate 
5i ~ 5 fi. X 10^^^ cm'^/s, g2 — 4: h x 10"^^ cm^/s, i.e. both nonlinear terms 
in Eq. ([T]) are positive (repulsive). However, in our case for the existence of self- 
trapped localized matter-waves we assume that the coefficient of cubic nonlinearity 
gi = Airh^as/m is shifted to negative (attractive) domain via change of the s-wave 
scattering length a^. by a Feshbach resonance technique. The strength of transverse 
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confinement used in matter-wave soliton experiments ^wj^ = 700-;- 800 provides 
a_L ~ 1 iim. Then, for the number of atoms in Rb BEC N ~ 10"* we get the fol- 
lowing values for the dimensionless coefficients of cubic and quintic nonlinear terms 
a = -giN/{2T:alhLu±) = 2(|as|/aj_)iV ~ 100, /3 = q2N^ /{^Tr'^a^^hujjL) - 200. In 
subsequent numerical simulations we employ the above estimates for the coefficients. 

Exact one soliton solutions of Eq. ((3]) in the absence of external potential 
{U{x) = 0) were found in Ref. For the case of self-focusing cubic {a > 0) 
and defocusing quintic (/3 > 0) nonlinearities, under normalization condition ([2]), 
the solution is 

/3a tanh(77) exp[z(ga; - /it)] Pip 1 rj 

V4/3 ^l + sech(r/)cosh(a:/a) V 3 ' a tanh(7y) ' 

where q and fi stand for the wave vector and chemical potential. In presence of a 
trap potential the Eq. ([3]) does not have analytic solution, and therefore one has to 
recourse to approximate methods. 

It is appropriate to mention, that cubic-quintic nonlinear Schrodinger equation 
([3| wi th flat-to p solutions appears in a variety of physical contexts in nonlinear 
optics ^ l ^^ l ^^ l, fluid dynamics -IZ,, plasma physics and BEC . Recently it was 
considered as a model equation describing bright solitons in the Tonks-Girargeau 
gas with dipolar interactions Therefore, the variational approach developed in 
this paper is of general interest for the above mentioned fields. 

3. Static variational approximation 

The variational approximation represents one of th e im portant theoretical tools 
for investigation of solitons in non-integrable models The success of VA essen- 
tially depends on the proper choice of a trial function. In particular, significant 
progress has been made with application of VA to nonlinear Schrodinger (NLS) 
type equations with cubic nonlinearity, employing Gaussian and hyperbolic secant 
trial functions. However, when the higher order nonlinear terms are included in 
the NLS equation, the shapes of localized states may significantly deviate from the 
above mentioned functions, and one has to consider other options. The possibility 
to perform analytic calculations is the major issue in selection of a trial function. 

For the NLS with competing cubic and quin tic nonlinearities , wh en the soliton 
features a fiat-top shape, a super-Gaussian and super-secant trial functions 
were shown to be appropriate for the description of self-trapping of laser beams in 
two-dimensional (2D) cubic-quintic nonlinear media. The behavior of soliton solu- 
tion of the NL S eq uation with arbitrary nonlinearity near the blow-up point was 
investigated in 1231 by means of VA based on a super-Gaussian trial function, and 
accurate estimate for the critical blow-up mass was found. Among other success- 
ful applications of the super-Gaussian trial function, description of the pulsating 
localized solutions of the cubic-quintic complex Ginzburg-Landau equation 1251 and 



stationary solutions to the NLS equation in a parabolic-index fiber can be men- 
tioned. 
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Although the Eq. ([3]) without external potential {U{x) — 0) has exact one soli- 
ton solution (|4]) , in Ref. [27| a variational approximation using the Gaussian trial 
function was developed for this case. Justification for the construction of VA when 
the exact soliton solution is available, can be found in some advantages provided by 
the VA in the analysis of existence and stability of solitons. For instance, when both 
of the nonlinear terms in Eq. ([3]) are focusing, the localized wave undergoes collapse 
(unlimited shrinking) if the norm of the wave function exceeds some critical value. 
The VA with a Gaussian trial function accurately predicts the corr esp onding thresh- 
old norm (the solution ceases to exist at this value of the norm) Furthermore, 
simple analytic relations between parameters of the localized state provided by VA 
facilitates its stability analysis by means of the Vakhitov - Kolokolov criterion 

However, the Gaussian trial function restricts the validity of the developed VA 
to specific areas of the parameter space of Eq. ([3]), since it is adequate only if both 
nonlinear terms are focusing, or when the effect of repulsive quintic term is weak 
compared to the attractive cubic one. In the opposite situation the localized solution 
acquires a flat-top shape, and one has to consider a different trial function. With this 
motivation in mind below we develop the VA for Eq. ([3]) using a super-Gaussian trial 
function, and apply it for the analysis of static and dynamic properties of flat-top 
matter- wave solitons. 



3.1. Stationary wave profile in a free space 

It is instructive to start with considering the variational solution of Eq. ([3]) in a 
free space and comparing the obtained wave profile with the available exact one 
soliton solution ([4]). The aim here is twofold. Form one side we can thereby check 
the accuracy of the VA, and from the other side, useful relations between soliton 
parameters will be obtained. 

In absence of external trap potential the governing equation has the following 
form 

#t + ^V'.. + a|V'l''/'-/3|V'lV = 0. (5) 

Although the coefficient of cubic nonlinearity a can be rescaled to one by trans- 
formations ip \fmli and /3 P/o? , we retain it for convenience of the time 
dependent VA to be considered later. The stationary states, which are looked for 
as ip{x,t) = 4>{x) exp(— i/it), satisfy the following equation 

M'/' + ^0.. + a'/'' - 13^' = 0. (6) 
The Lagrangian density generating this equation is 

Since the typical localized solutions of Eq. ([5]) for competing nonlinearities (at- 
tractive cubic and repulsive quintic) are the "flat-top" solitons, we employ a super- 
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Gaussian trial function 



(x) — A exp 



1 / 



2 \aJ 



(8) 



with A, a, m being the variational parameters, corresponding to the amplitude, 
width and super-Gaussian index of the soliton, respectively. 

The averaged Lagrangian L = Cdx computed with this ansatz is 



L 



r(2 -M)-[ A^at, + - ) r(l + M), 



8aAf 



(9) 



where r(a;) is the gamma function, and the notation M ~ l/2m is introduced. 
The stationary values for the amplitude {A), width (a) and reduced super-Gaussian 
index [M) of the wave profile are found from the variational equations 



dL _ dL _ dL 
9^2 Ih.^^' dM 



0. 



(10) 



Straightforward calculations, with taking into regard the above mentioned normal- 
ization 

/ + 00 
(j)^{x)dx^2A'^ar{l + M) = l, (11) 
-CO 

yield the following relations between variational parameters 



(3/2) 



M+l 



2/3-^3^^+1 r(M) r(2- Af) 
1 



1/2 



2A2r(i + M) ' 

2,aA^ 2(3A'^ 

2M+2 3M+T' 



(12) 
(13) 
(14) 



r(M)r(2 - M) [M-i + 1^(2 - M) - + M) - 21n2] + /3^|^^ 



0, (15) 



where il^{x) = ^lnr(a:) is the digamma function'^ (not to be confused with the 
wave function). The Eqs. p2 |l - (fT5|) are sufficient to determine the four parameters 
of the localized state (A, a, to, /j,). An example of stationary solution of Eq. ([5]) for 
a particular set of parameters is presented in Fig. [T] As can be seen from this figure 
the agreement between the exact and VA wave profiles for the flat-top soliton is 
quite good. 

An important observation following from the above analysis is that, in the flat- 
top regime the super-Gaussian index M (and therefore m = 1/2^/) for the soliton 
does not depend on the coefficient of cubic nonlinearity a. In fact this is the mani- 
festation of above mentioned rescaling property of Eq. ([5]) . This property later will 
be used in construction of the time-dependent VA. In the next subsection we extend 
the static VA for the case when the fiat-top soliton is confined to a harmonic trap. 
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Fig. 1. (Color online) Comparison of the exact solution (red continuous line) given by Eq. Q 
with the prediction of VA (blue dashed line) for the flat-top soliton according to Eqs. I)12|l - I|15| l 
for a = 100, /3 = 200. Parameters of the flat-top soliton found from VA are A = 0.62, a = 1.41, 
M = 0.155, At = -9.35, m = 1/2M = 3.22. 



3.2. Flat-top soliton in harmonic trap 

Usual experimental settings for BEC involve magnetic or optical traps designed for 
holding and manipulation with matter-waves. In a confining trap potential the soli- 
ton experiences deformation to some extent with respect to the free space condition. 
Below we consider the harmonic trap potential U(x) = ex^ in Eq. ([3]) and derive 
corresponding VA equations. 

The Lagrangian density generating the Eq. ([3]) is 



2^ -2^(--)^^-I^ ' 6- 



(16) 



and the corresponding averaged Lagrangian computed with the super-Gaussian 
ansatz ^ has the form 



L : 



Ah 



A'^af] 



T)r(i+M). (17) 



8aM 3 V" ' 2^^+i 3J^^^ 

Application of the conditions (fTU]) to this averaged Lagrangian yields the following 
relations between parameters of the flat-top soliton 

2a2er(l + 3A/) 3aA^ 213 A'^ 



3r(l + M) 2^^+2 3''^+^' 
The effective width of the soliton is found as the root of the algebraic equation 



(18) 



a + pa + q = 0, 



(19) 



with the coefficients 



3a 



2A/+3r(H-3M)£' 



2/3 + 3*=f+ir(M)r(2 - M) 
8 • 3A^r(l + M)r(l + 3M) e ' 
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Real and positive root of the Eq. is 



I 2p ^/i 



(20) 



where 



6r 




The counterpart of the Eq. (fT5|) for the trapped sohton has the form 

e r(l + MI) (2V'(1 + M) - Stpil + 3M)) a a {2 ln2 + V(l + M)) 



3 2*f+3 
r(2 - M) (1 + MV'(2 - M)) 13 (hi3 + ipil + M)) 



0. (21) 



8M2 4 • 3A^+ir(l + M) 

In order to find the shape of the flat-top sohton conflned to a parabohc trap at flrst 
we solve this equation with respect to M, substituting a from Eq. (|20p . Subsequently 
the width a is computed from Eq. ([^0]) using the value of M found as a root of 
Eq. (|2ip . Next, the amplitude A is computed from the expression for the norm 
(fTTj) . and chemical potential is found from Eq. In the left panel of Fig. [5] we 

illustrate the shapes of flat-top solitons for two strengths of the trap potential e 
as predicted by VA, and as found from the original GPE ([3]) by imaginary time 
relaxation method . As can be seen from this figure, stronger parabolic trap leads 




to more deformation of the soliton shape. Specifically, under the effect of a parabolic 
trap the flat-top soliton shrinks, while its amplitude increases. Dependence of the 
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width (a) and super- Gaussian index (m) as a function of the strength of the trap 
potential (e) is depicted in the right panel of Fig. [3l The deviation between the 
results of numerical solution of the original GPE ([3]) in imaginary time 1^ and 
prediction of the VA increases as the trap potential becomes stronger. However, the 
discrepancy for the given range of e is less than 3 %, thus the agreement is quite 
good. 




Fig. 3. (Color online) Left panel: The width of the flat-top soliton as a function of the strength of 
the trap potential. Red solid line corresponds to numerical solution of the original GPE ||3]l, while 
blue dashed line is the prediction of VA for a = 100 and fS = 400. Right panel: The super-Gaussian 
index in Eq. ^ as a function of the strength of the trap potential. 



4. Time dependent variational approximation 

The time dependent VA allows to investigate soliton dynamics under external or 
internal perturbations, such as varying strength of the trap potential, chirp imprint- 
ing or alternating coefficient of nonlinear interaction. Below we develop the dynamic 
version of VA for the latter case, when the coefficient of cubic nonlinearity in Eq. 
([1]) is periodically varied in time, via the s-wave scattering length, by means of a 
Feshbach resonance technique I^. This setting is frequently called as the nonlinearity 
management 

For simplicity we consider the case when the trap potential is absent. In dimen- 
sionless units the governing equation has the form 

#t + ^^cc. + aWIV-lV - /^IV'lV = 0, (22) 

where a{t) — —2as{t)N/a± is the time dependent strength of the cubic nonlinearity. 
In the following we consider the periodic a{t) — ao[l + dsm{ujt)] variation of this 
parameter, with the magnitude 6 and frequency around its stationary value uq. 
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To construct the time dependent VA we consider the fohowing Lagrangian den- 
sity, generating the Eq. (|22|) 



c = '-m: - r^t) + - ^iv'i' + f t)f. (23) 

While selecting the trial function we recall that for the flat-top soliton the super- 
Gaussian index m does not explicitly depend on the strength of cubic nonlinearity 
a, according to Eq. Therefore, when the parameter a is subject to variations, 
the following trail function with constant to is appropriate 

2m 



t) = A{£) exp 



'2 \j4tj 



(24) 



where A(t),a{t),b{t) and 0(i) are the time dependent variational parameters, de- 
noting the amplitude, width, chirp and phase of the flat-top soliton, respectively. 

Substituting the ansatz into Eq. (|23p and performing spatial integration 
yields the following averaged Lagrangian 



2,2^r(i + 3M) , r(2-M) 



3r(l-hM) 8Afr(l + M), 

a{t) 13 



2*^+2r(l + M)a ' 4-3A^+ir2(l + M)a2' ^^^^ 

where we have used the normalization condition (jlip to eliminate A, and employed 
the notation M — 1/2to. Variation of this Lagrangian with respect to time depen- 
dent parameters a, h and (j) leads to the following equation for the width of the 
flat-top soliton 

h{M)+ph{M) a{t)h{M) 
an = n 5 , (^Dj 



where 



MM) 



4Mr(l + 3Af)' 2-3^^r(l + M)r(l-h3M)' 

3 



2A/+2r(i + 3Af)' 

This equation is analogous to the equation of motion of a unit mass particle in 
anharmonic potential att = —dU{a)/da, with 

^ ^ MM) + (3 MM) _ aoMM) ^ 

where ao is the stationary value of the cubic nonlinear coefflcient. The shape of 
this potential is depicted in the left panel of Fig. 21 The minimum of the poten- 
tial corresponds to the width of the unperturbed flat-top soliton oq, found from 
Eqs. ([T^ - P3)) . Weakly perturbed (chirped) soliton performs small amplitude oscil- 
lations around the minimum of the potential (I27p with a frequency 

c^o = [3(/i(M) + PMM)) - 2aoah{M)f/^/al. (28) 
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Oscillations of the mean square width of a flat-top sohton, defined as 




Fig. 4. (Color online) Left panel: The potential 1271 1 for the effective particle model Eq. II26I I 
for two sets of parameters found from stationary version of VA: a = 100, /3 = 200, M = 0.155, 
ao = 1.41 and cx = 100, /3 = 400, M = 0.111, ao = 2.78. The minimum of the potential corresponds 
to the stationary width (ao) of the flat-top soliton. Right panel: Oscillation of the amplitude under 
slowly varying strength of cubic nonlinearity a{t) = ao[l + 5 sm{wt)\, with 5 = 0.05, u) = 0.1, found 
from numerical solution of the GPE l|22| l (red solid line), and as predicted by VA Eq. I|26| l (blue 
dashed line). 

oo 

-^(i)= / x^\i.{xMdx^a\t)^^±^^, (29) 

— oo 

is shown in the right panel of Fig. |4l Apparently, the dynamics described by the VA 
Eq. (|26|) agrees quite well with numerical solution of the GPE (p2)) . The peculiarity 
of the dynamics is that, the fine structure due to fast internal vibrations of the 
soliton with a frequency (j28p is superimposed upon the slow dynamics under non- 
linearity management. For the parameter settings of the right panel of Fig. 21 the 
frequency of internal vibrations found from numerical solution of the GPE (j22|) is 
070 — 1.76, while the prediction of the VA according to Eq. is loq = 1.89. When 
the frequency of nonlinearity management is close to this frequency, the dynamics is 
highly irregular. Surface plot of the flat-top soliton evolving under slow nonlinearity 
management is shown in Fig. [5l 

Numerical simulations of the GPE ([22| are performed by means of the split-step 
fast-Fourier-transform method ' '^ ' in a spatial domain of length L = Stt with 
2048 modes. The time step was At = 0.001. To control the numerical results, we 
monitored the accuracy of normalization condition ([2]), which showed conservation 
during the time evolution to precision better than 10~^ in normalized units. To 
prevent re-entering of the linear waves emitted by the perturbed soliton into the 
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integration domain, absorbers were installed at domain boundaries. To obtain the 
stationary profiles of flat-top solitons in external potentials for Eq. ([3]) we employed 
the imaginary-time relaxation method for finding ground states in NLSE-based 

I on I 

models 




Fig. 5. (Color online) Surface plot of a flat-top soliton evolving under nonlinearity management 
according to numerical solution of GPE II22I I. Parameter settings correspond to the right panel of 
Fig.E] 



5. Conclusions 

We have developed a variational approach for flat-top solitons, and employed it 
for the analysis of static and dynamic properties of matter- wave solitons which can 
exist in EEC's with significant contribution of three-body interactions. The accuracy 
of the developed approach has been verified by comparing the predictions of VA 
equations with corresponding data from the numerical solution of the governing 
GPE, and good agreement is found between the two results. Although the emphasis 
was given to BEC applications, the developed theory is general and may apply to 
nonlinear optics phenomena in materials with cubic-quintic nonlinearity. 
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